# Replikationsdaten f�r
# Weidmann, Nils B. (2009)
# "Kausalprozesse auf mehreren Ebenen:
# Ein Pl�doyer f�r vertikale Integration
# in der B�rgerkriegsforschung"
# Politische Vierteljahresschrift, 2009

# Diese Analyse wurde mit R Version 2.6.2
# und dem lme4 package, Version 0.99875-9
# durchgef�hrt.

library(lme4)

# Datensatz einlesen
units <- read.table("pvs.txt", header=T)

# Hierarchisches Modell
model2 <- lmer(as.factor(CONFLICT) ~ ethfrac + gdpen + polity2 + elf + diamonds + terrain + log10(population) + W.CONFLICT + (1 |ccode),  data=units, family=binomial(link="probit"))
summary(model2)

# eigener data.frame f�r L�nder
cdata <- subset(data.frame(units), select=c(ccode, ethfrac, polity2, gdpen, EVCNT, CNTRY_NAME))
for (c in unique(cdata$ccode)) {
  confcnt <- max(cdata[cdata$ccode==c,]$EVCNT)
  cdata$conflict[cdata$ccode==c] <- ifelse(confcnt>0,1,0)
}
cdata <- cdata[!duplicated(cdata$ccode),]
cdata <- cdata[order(cdata$ccode),]

# Abb. 4 links
gdp.mean <- mean(units$gdpen)
polity2.mean <- mean(units$polity2)
a.hat.model2 <- fixef(model2)[1] + fixef(model2)[2]*cdata$ethfrac + fixef(model2)[3]*gdp.mean + fixef(model2)[4]*polity2.mean + as.vector(ranef(model2)$ccode)
plot(cdata$ethfrac, t(a.hat.model2), cex=2, bty="n", xlim=c(0, 1), ylim=c(-20, -12), xlab="Ethnische Fragmentierung (auf nationaler Ebene)", ylab=expression(paste("Gesch�tzes Intercept","  ", alpha[j])))
abline(fixef(model2)[1] + fixef(model2)[3]*gdp.mean + fixef(model2)[4]*polity2.mean , fixef(model2)[2])
for (i in 1:nrow(cdata)) {
  if (cdata[i,]$conflict) {
     points(cdata$ethfrac[i], t(a.hat.model2)[i], pch=19, cex=2)
  }
}

# Abb. 4 rechts
invlogit <- function(x) {
    exp(x) / (1+exp(x))
}

plot(1, type="n", xlim=c(0,1), ylim=c(0.4,1), bty="n", xlab="Ethnische Fragmentierung (auf regionaler Ebene)", ylab="P(Konflikt)")
f <- fixef(model2)
linetypes <- c(1,2,3)
ind <- 0
for (c in c(2,3,6)) {
  ind <- ind+1
  ccode <- cdata[c,]$ccode
  diamonds.mean <- mean(units$diamonds[units$ccode==ccode])
  terrain.mean <- mean(units$terrain[units$ccode==ccode])
  pop.mean <- mean(log10(units$population[units$ccode==ccode]))
  W.conflict.mean <- mean(units$W.CONFLICT[units$ccode==ccode])
  r <- cdata[c,]
  curve(invlogit(f[1] + f[2]*r$ethfrac + f[3]*r$gdpen + f[4]*r$polity2 + data.frame(ranef(model2)[1])[c,] + f[5]*x + f[6]*diamonds.mean + f[7]*terrain.mean + f[8]*pop.mean + f[9]*W.conflict.mean), add=T, lty=linetypes[ind], lwd=2)
}
text(0.4, 0.92, "Kongo (DRC)")
text(0.4, 0.55, "Uganda")
text(0.4, 0.72, "Rwanda")


